#
#The Following R-Code runs OP, ZiOP, and ZiOPC models of dyadic hosility,
#using datasets similar to those used by Senese (1997).  
#
#These correspond to Table 2 in our paper
#
#
#All OP, ZiOP, and ZiOPC models are then compared with AIC statistics, 
#Vuong tests, and likelihood ratio tests (these are referenced in the text of our paper).
#
#
#

#clear memory
rm( list=ls() ) 

## set the memory to the max
memory.limit(size=4095)
memory.size(max = TRUE)
						      
#install.packages
library(mvtnorm)
library(car)                                                        
library(scatterplot3d)
library(plotrix)
library(Hmisc)
library(Zelig)
library(foreign, pos=4)
library(MASS)
library(maxLik)

#set the seed
set.seed(2)

#read in formatted Senese replication dataset
senese<- read.dta("replication.dta", 
  convert.dates=TRUE, convert.factors=TRUE, missing.type=TRUE, 
  convert.underscore=TRUE, warn.missing.labels=TRUE)

#summarize dataset
summary(senese)

#subset the dataset to include only necessary variables
newdata.onset<-senese[c(1,2,3,4,5,6,7,8,9,10,11,12,13)]

#omit observations with missing data from dataset
newdata.onset <- na.omit(newdata.onset)

#summarize new dataset to ensure missignness has been removed
summary(newdata.onset)

#create a second datset with no peace-year observations, as used in Senese (1997)
nz.newdata.onset <- subset(newdata.onset, newdata.onset$s.hostility!=0)



#Begin Analysis#

############################################################################
########################Truncated OP Model##################################
############################################################################

#This first model estimates an OP model of hostility when all (zero) peace-years
#are removed, which is consistent with Senese (1997).

#write a function to estimate an OP model 
OP <- function(b,data) {					      
  y<-data[,13]			# hostility
  x1<-data[,6]			# contig
  x2<-data[,7]			# alliance
  x3<-data[,8]			# dem
  x4<-data[,9]			# power		                                             
  x5<-data[,11]			# maturity one
  x6<-data[,12]			# maturity none                                               	                                             					     	      
  z<-cbind(x1,x2,x3,x5,x6)	# covariates 
  tau<- rep(0,6) 		                              
  tau[1]<--(Inf)
  tau[2]<- b[1]
  tau[3]<- tau[2]+exp(b[2])
  tau[4]<- tau[3]+exp(b[3])
  tau[6]<- (Inf)
  n=nrow(data)							      							      
  llik <- matrix(0, nrow=n, ncol = 1)   			      					      					              					      				      
  ZG  <- b[4]*z[,1] + b[5]*z[,2] + b[6]*z[,3] + b[7]*z[,4]+ b[8]*z[,5]		              
  llik<-ifelse(y==1,log((pnorm(tau[2]-ZG))),
	ifelse(y==2,log((pnorm(tau[3]-ZG) - pnorm(tau[2]-ZG))),
        ifelse(y==3,log((pnorm(tau[4]-ZG) - pnorm(tau[3]-ZG))),
        log((1 - pnorm(tau[4]-ZG)))     
	)))		
        llik<--1*sum(llik)					      
    	return(llik)
}										      

#give the OP model resonable starting parameters
b<-rbind(2.899,-4.34693, -2.49664, 1.21145,-0.12,-0.44,0.54,0.02)

#optimize the OP model
output.OP<-nlm(f=OP, p=b, data=nz.newdata.onset, iterlim = 500, hessian=TRUE)

#show the output to get the coefficient estimates
print(output.OP)

#recover the variance covariance matrix
op.vcv<-solve(output.OP$hessian)

#show the variance covariance matrix to recover the standard errors
op.vcv



############################################################################
########################Non-Truncated OP Model##############################
############################################################################

#This model estimates an OP model of hostility when peace-years are included.

#write a function to estimate an OP model 
OP <- function(b,data) {					      
  y<-data[,13]			# hostility
  x1<-data[,6]			# contig
  x2<-data[,7]			# alliance
  x3<-data[,8]			# dem
  x4<-data[,9]			# power		                                             
  x5<-data[,11]			# maturity one
  x6<-data[,12]			# maturity none                                                     	                                             					     	      
  z<-cbind(x1,x2,x3,x5,x6)	# covariates 
  tau<- rep(0,6) 		                              
  tau[1]<--(Inf)
  tau[2]<- b[1]
  tau[3]<- tau[2]+exp(b[2])
  tau[4]<- tau[3]+exp(b[3])
  tau[5]<- tau[4]+exp(b[4])
  tau[6]<- (Inf)
  n=nrow(data)							      							      
  llik <- matrix(0, nrow=n, ncol = 1)   			      					      					              					      				      
  ZG  <- b[5]*z[,1] + b[6]*z[,2] + b[7]*z[,3] + b[8]*z[,4] + b[9]*z[,5]	
  llik<-ifelse(y==0,log((pnorm(tau[2]-ZG))),
	ifelse(y==1,log((pnorm(tau[3]-ZG) - pnorm(tau[2]-ZG))),
	ifelse(y==2,log((pnorm(tau[4]-ZG) - pnorm(tau[3]-ZG))),
	ifelse(y==3,log((pnorm(tau[5]-ZG) - pnorm(tau[4]-ZG))),
	log((1 - pnorm(tau[5]-ZG)))
	))))

        llik<--1*sum(llik)					      
    	return(llik)
}								      

#give the OP model resonable starting parameters
b<-rbind(2.899,-4.34693, -2.49664, -0.54123, 1.21145,-0.12,-0.44,0.54,0.02)

#optimize the OP model
output.OP<-nlm(f=OP, p=b, data=newdata.onset, iterlim = 500, hessian=TRUE)

#show the output to get the coefficient estimates
print(output.OP)

#recover the variance covariance matrix
op.vcv<-solve(output.OP$hessian)

#show the variance covariance matrix to recover the standard errors
op.vcv




############################################################################
############################Main ZiOP Model#################################
############################################################################

#This model estimates our main ZiOP model of hostility when peace-years are included.

#write a function to estimate a ZiOP model 
ZIOP <- function(b,data) {					      
  y<-data[,13]			# hostility
  x1<-data[,6]			# contig
  x2<-data[,7]			# alliance
  x3<-data[,8]			# dem
  x4<-data[,9]			# power		                                             
  x5<-data[,11]			# maturity one
  x6<-data[,12]			# maturity none
  x<-cbind(1,x1,x2,x3,x4)	# inflation stage covariates                                                   	                                             					     	      
  z<-cbind(x1,x2,x3,x5,x6)	# outcome stage covariates  

  tau<- rep(0,6) 		                              
  tau[1]<--(Inf)
  tau[2]<- b[1]
  tau[3]<- tau[2]+exp(b[2])
  tau[4]<- tau[3]+exp(b[3])
  tau[5]<- tau[4]+exp(b[4])
  tau[6]<- (Inf)
  n=nrow(data)							      							      
  llik <- matrix(0, nrow=n, ncol = 1)   			      					      					              					      						      
  XB  <-b[5]*x[,1] + b[6]*x[,2] + b[7]*x[,3] +b[8]*x[,4]+b[9]*x[,5]					      
  ZG  <- b[10]*z[,1] + b[11]*z[,2] +b[12]*z[,3] +b[13]*z[,4] +b[14]*z[,5]  	

  llik<-ifelse(y==0,log( (1-pnorm(XB))+ (pnorm(XB)) * ( pnorm(tau[2]-ZG)) ),
	ifelse(y==1,log( (pnorm(XB)) * (pnorm(tau[3]-ZG) - pnorm(tau[2]-ZG))),
	ifelse(y==2,log( (pnorm(XB)) * (pnorm(tau[4]-ZG) - pnorm(tau[3]-ZG))),
	ifelse(y==3,log( (pnorm(XB)) * (pnorm(tau[5]-ZG) - pnorm(tau[4]-ZG))),
	log( (pnorm(XB)) * (1-pnorm(tau[5]-ZG)))
	))))

        llik<--1*sum(llik)					      
    	return(llik)
}


#give the ZiOP model resonable starting parameters
b<-rbind(2.899,-4.34693, -2.49664, -0.54123,-2.1,1.21,-0.12,2.44,.1,0.54,0.02,1.21,-0.12,.1)

#optimize the ZiOP model
output.ZIOP<-nlm(f=ZIOP, p=b, iterlim = 500, data=newdata.onset, hessian=TRUE)

#show the output to get the coefficient estimates
print(output.ZIOP)

#recover the variance covariance matrix
output.ZIOP.vcv<-solve(output.ZIOP$hessian)

#show the variance covariance matrix to recover the standard errors
output.ZIOP.vcv




############################################################################
########################ZIOP Robustness Model###############################
############################################################################

#This model estimates a ZiOP model of hostility when peace-years are included, and when
#both maturity variables are additionally included in the inflation stage.

#write a function to estimate a ZiOP model
ZIOP <- function(b,data) {					      
  y<-data[,13]			# hostility
  x1<-data[,6]			# contig
  x2<-data[,7]			# alliance
  x3<-data[,8]			# dem
  x4<-data[,9]			# power		                                             
  x5<-data[,11]			# maturity one
  x6<-data[,12]			# maturity none
  x<-cbind(1,x1,x2,x3,x4,x5,x6)	# inflation stage covariates                                                    	                                             					     	      
  z<-cbind(x1,x2,x3,x5,x6)	# outcome stage covariates  
  tau<- rep(0,6) 		                              
  tau[1]<--(Inf)
  tau[2]<- b[1]
  tau[3]<- tau[2]+exp(b[2])
  tau[4]<- tau[3]+exp(b[3])
  tau[5]<- tau[4]+exp(b[4])
  tau[6]<- (Inf)
  n=nrow(data)							      							      
  llik <- matrix(0, nrow=n, ncol = 1)   			      					      
  XB  <-b[5]*x[,1] + b[6]*x[,2] + b[7]*x[,3] +b[8]*x[,4]+b[9]*x[,5]+b[10]*x[,6]+b[11]*x[,7]						      
  ZG  <- b[12]*z[,1] + b[13]*z[,2] +b[14]*z[,3] +b[15]*z[,4] + b[16]*z[,5]				              					      
	              

  llik<-ifelse(y==0,log( (1-pnorm(XB))+ (pnorm(XB)) * ( pnorm(tau[2]-ZG)) ),
	ifelse(y==1,log( (pnorm(XB)) * (pnorm(tau[3]-ZG) - pnorm(tau[2]-ZG))),
	ifelse(y==2,log( (pnorm(XB)) * (pnorm(tau[4]-ZG) - pnorm(tau[3]-ZG))),
	ifelse(y==3,log( (pnorm(XB)) * (pnorm(tau[5]-ZG) - pnorm(tau[4]-ZG))),
	log( (pnorm(XB)) * (1-pnorm(tau[5]-ZG)))
	))))

        llik<--1*sum(llik)					      
    	return(llik)
}


#give the ZiOP model resonable starting parameters
b<-rbind(2.899,-4.34693, -2.49664, -0.54123,-2.1,1.21,-0.12,2.44,.1,0.54,0.02,1.21,-0.12,-0.44,-0.12,-0.44)

#optimize the ZiOP model
output.ZIOP<-nlm(f=ZIOP, p=b, iterlim = 500, data=newdata.onset, hessian=TRUE)

#show the output to get the coefficient estimates
print(output.ZIOP)

#recover the variance covariance matrix
output.ZIOP.vcv<-solve(output.ZIOP$hessian)

#show the variance covariance matrix to recover the standard errors
output.ZIOP.vcv




#############################################################################
###########################Main ZIOPC Model##################################
#############################################################################

#This model estimates our main ZiOPC model of hostility when peace-years are included.

#write a function to estimate a ZiOPC model
ZIOPC <- function(b,data) {					      
  y<-data[,13]			# hostility
  x1<-data[,6]			# contig
  x2<-data[,7]			# alliance
  x3<-data[,8]			# dem
  x4<-data[,9]			# power		                                             
  x5<-data[,11]			# maturity one
  x6<-data[,12]			# maturity none
  x<-cbind(1,x1,x2,x3,x4)	# inflation stage covariates   	                                                 	                                             					     	      
  z<-cbind(x1,x2,x3,x5,x6)	# outcome stage covariates  	
  tau<- rep(0,6) 		                              
  tau[1]<--(Inf)
  tau[2]<- b[1]
  tau[3]<- tau[2]+exp(b[2])
  tau[4]<- tau[3]+exp(b[3])
  tau[5]<- tau[4]+exp(b[4])
  tau[6]<- (Inf)
  n=nrow(data)								      							      
  llik <- matrix(0, nrow=n, ncol = 1)   			      					      
  for(i in 1:n){					              					      
        	B<-b[5:9]						      
 		XB  <- B %*% x[i,] 				      
        	G<-b[10:14]
       		ZG  <- G %*% z[i,]   
                rho<-b[15] 
                sigma<-matrix(c(1,rho,rho,1), nrow=2, ncol=2)
                negsigma<-matrix(c(1,-rho,-rho,1), nrow=2, ncol=2)		              
                if(identical(y[i],0)){llik[i]<-log((1-pnorm(XB))+ pmvnorm(lower=c(-Inf,-Inf), upper=c(XB,(tau[2]-ZG)),mean=c(0,0),sigma=negsigma))}
                else if(identical(y[i],1)){llik[i]<-log(pmvnorm(lower=c(-Inf,-Inf), upper=c(XB,(tau[3]-ZG)),mean=c(0,0),sigma=negsigma)-pmvnorm(lower=c(-Inf,-Inf), upper=c(XB,(tau[2]-ZG)),mean=c(0,0),sigma=negsigma))}
		else if(identical(y[i],2)){llik[i]<-log(pmvnorm(lower=c(-Inf,-Inf), upper=c(XB,(tau[4]-ZG)),mean=c(0,0),sigma=negsigma)-pmvnorm(lower=c(-Inf,-Inf), upper=c(XB,(tau[3]-ZG)),mean=c(0,0),sigma=negsigma))}
		else if(identical(y[i],3)){llik[i]<-log(pmvnorm(lower=c(-Inf,-Inf), upper=c(XB,(tau[5]-ZG)),mean=c(0,0),sigma=negsigma)-pmvnorm(lower=c(-Inf,-Inf), upper=c(XB,(tau[4]-ZG)),mean=c(0,0),sigma=negsigma))}
	        else if(identical(y[i],4)){llik[i]<-log(pmvnorm(lower=c(-Inf,-Inf), upper=c(XB,(ZG-tau[5])),mean=c(0,0),sigma=sigma))}     
		}
        llik<--1*sum(llik)					      
    	return(llik)
}								      

############################################################################

#give the ZiOPC model resonable starting parameters
b<-rbind(0.402, -3.8, -2, -0.22, -2.63, 3.17, 0.41, -0.12, 0.30, -1.02, -0.26, -0.44, -0.06, .1,-0.122)

#optimize the ZiOPC model
output.ZIOPC<-nlm(f=ZIOPC, p=b, iterlim = 500, data=newdata.onset, hessian=TRUE)

#show the output to get the coefficient estimates
print(output.ZIOPC)

#recover the variance covariance matrix
ziopc.vcv<-solve(output.ZIOPC$hessian)

#show the variance covariance matrix to recover the standard errors
ziopc.vcv 


############################################################################
########################ZIOPC Robustness Model##############################
############################################################################

#This model estimates a ZiOPC model of hostility when peace-years are included, and when
#both maturity variables are additionally included in the inflation stage.

#write a function to estimate a ZiOPC model
ZIOPC <- function(b,data) {					      
  y<-data[,13]			# hostility
  x1<-data[,6]			# contig
  x2<-data[,7]			# alliance
  x3<-data[,8]			# dem
  x4<-data[,9]			# power		                                             
  x5<-data[,11]			# maturity one
  x6<-data[,12]			# maturity none
  x<-cbind(1,x1,x2,x3,x4,x5,x6)	# inflation stage covariates                                                    	                                             					     	      
  z<-cbind(x1,x2,x3,x5,x6) 	# outcome stage covariates 
  tau<- rep(0,6) 		                              
  tau[1]<--(Inf)
  tau[2]<- b[1]
  tau[3]<- tau[2]+exp(b[2])
  tau[4]<- tau[3]+exp(b[3])
  tau[5]<- tau[4]+exp(b[4])
  tau[6]<- (Inf)
  n=nrow(data)								      							      
  llik <- matrix(0, nrow=n, ncol = 1)   			      					      
  for(i in 1:n){					              					      
        	B<-b[5:11]						      
 		XB  <- B %*% x[i,] 				      
        	G<-b[12:16]
       		ZG  <- G %*% z[i,]   
                rho<-b[17] 
                sigma<-matrix(c(1,rho,rho,1), nrow=2, ncol=2)
                negsigma<-matrix(c(1,-rho,-rho,1), nrow=2, ncol=2)		              
                if(identical(y[i],0)){llik[i]<-log((1-pnorm(XB))+ pmvnorm(lower=c(-Inf,-Inf), upper=c(XB,(tau[2]-ZG)),mean=c(0,0),sigma=negsigma))}
                else if(identical(y[i],1)){llik[i]<-log(pmvnorm(lower=c(-Inf,-Inf), upper=c(XB,(tau[3]-ZG)),mean=c(0,0),sigma=negsigma)-pmvnorm(lower=c(-Inf,-Inf), upper=c(XB,(tau[2]-ZG)),mean=c(0,0),sigma=negsigma))}
		else if(identical(y[i],2)){llik[i]<-log(pmvnorm(lower=c(-Inf,-Inf), upper=c(XB,(tau[4]-ZG)),mean=c(0,0),sigma=negsigma)-pmvnorm(lower=c(-Inf,-Inf), upper=c(XB,(tau[3]-ZG)),mean=c(0,0),sigma=negsigma))}
		else if(identical(y[i],3)){llik[i]<-log(pmvnorm(lower=c(-Inf,-Inf), upper=c(XB,(tau[5]-ZG)),mean=c(0,0),sigma=negsigma)-pmvnorm(lower=c(-Inf,-Inf), upper=c(XB,(tau[4]-ZG)),mean=c(0,0),sigma=negsigma))}
	        else if(identical(y[i],4)){llik[i]<-log(pmvnorm(lower=c(-Inf,-Inf), upper=c(XB,(ZG-tau[5])),mean=c(0,0),sigma=sigma))}     
		}
        llik<--1*sum(llik)					      
    	return(llik)
}								      

############################################################################

#give the ZiOPC model resonable starting parameters
b<-rbind(0.402, -3.8, -2, -0.22, -2.63, 3.17, 0.41, -0.12, 0.30, 0.37,.1, -1.02, -0.26, -0.44, -0.06, -0.122, -.2)

#optimize the ZiOPC model
output.ZIOPC<-nlm(f=ZIOPC, p=b, iterlim = 500, data=newdata.onset, hessian=TRUE)

#show the output to get the coefficient estimates
print(output.ZIOPC)

#recover the variance covariance matrix
ziopc.vcv<-solve(output.ZIOPC$hessian)

#show the variance covariance matrix to recover the standard errors
ziopc.vcv 



############################################################################
#############################AIC Statistics#################################
############################################################################

#Truncated Ordered Probit
AIC.OP1<--2*(-1984)+2*(8)
AIC.OP1

#Non-Truncated Ordered Probit
AIC.OP2<--2*(-12773)+2*(9)
AIC.OP2

#Main ZiOP Model
AIC.ZIOP1<--2*(-11901)+2*(14)
AIC.ZIOP1

#ZiOP Robustness Model
AIC.ZIOP2<--2*(-11892)+2*(16)
AIC.ZIOP2

#Main ZiOPC Model
AIC.ZIOPC1<--2*(-11901)+2*(15)
AIC.ZIOPC1

#ZiOPC Robustness Model
AIC.ZIOPC2<--2*(-11891)+2*(17)
AIC.ZIOPC2


############################################################################
###############################Vuong Tests##################################
############################################################################

#First, compare main ZiOP(C) models to non-truncated OP model

#set the parameters
  y<-newdata.onset[,13]		# hostility
  x1<-newdata.onset[,6]		# contig
  x2<-newdata.onset[,7]		# alliance
  x3<-newdata.onset[,8]		# dem
  x4<-newdata.onset[,9]		# power		                                             
  x5<-newdata.onset[,11]	# maturity one
  x6<-newdata.onset[,12]	# maturity none
  x<-cbind(1,x1,x2,x3,x4)       # inflation stage covariates                                            	                                             					     	      
  z<-cbind(x1,x2,x3,x5,x6) 	# outcome stage covariates
  n<-nrow(newdata.onset)	# number of observations

     
#Vuong test comparing Main ZIOP to Non-Truncated OP
m<-matrix(0,nrow=n,ncol=1)
  for(i in 1:n){					              					      
	ZG.OP<-(1.18620*z[i,1]+-0.01618*z[i,2]+-0.32590*z[i,3]+0.01053*z[i,4]+-0.14820*z[i,5])
	cut1.OP<-2.72527
	cut2.OP<-2.72527+exp(-4.28509)
	cut3.OP<-2.72527+exp(-4.28509)+exp(-2.46453)
	cut4.OP<-2.72527+exp(-4.28509)+exp(-2.46453)+exp(-0.65400)

	XB.ZIOP<-(-2.73269*x[i,1]+3.74966*x[i,2]+0.28252*x[i,3]+-0.15223*x[i,4]+1.00665*x[i,5])
	ZG.ZIOP<-(-1.01742*z[i,1]+-0.21067*z[i,2]+-0.44175*z[i,3]+0.09286*z[i,4]+0.06755*z[i,5])
	cut1.ZIOP<-0.51352
	cut2.ZIOP<-0.51352+exp(-3.78640)
	cut3.ZIOP<-0.51352+exp(-3.78640)+exp(-1.98420)
	cut4.ZIOP<-0.51352+exp(-3.78640)+exp(-1.98420)+exp(-0.20561)
          
        if((y[i]==0)){m[i]<-log(pnorm(cut1.OP-ZG.OP)/((1-pnorm(XB.ZIOP))+(pnorm(XB.ZIOP))*(pnorm(cut1.ZIOP-ZG.ZIOP))))}
        else if((y[i]==1)){m[i]<-log((pnorm(cut2.OP-ZG.OP)-pnorm(cut1.OP-ZG.OP))/((pnorm(XB.ZIOP))*((pnorm(cut2.ZIOP-ZG.ZIOP)) - (pnorm(cut1.ZIOP-ZG.ZIOP)))))}
        else if((y[i]==2)){m[i]<-log((pnorm(cut3.OP-ZG.OP)-pnorm(cut2.OP-ZG.OP))/((pnorm(XB.ZIOP))*((pnorm(cut3.ZIOP-ZG.ZIOP)) - (pnorm(cut2.ZIOP-ZG.ZIOP)))))}
	else if((y[i]==3)){m[i]<-log((pnorm(cut4.OP-ZG.OP)-pnorm(cut3.OP-ZG.OP))/((pnorm(XB.ZIOP))*((pnorm(cut4.ZIOP-ZG.ZIOP)) - (pnorm(cut3.ZIOP-ZG.ZIOP)))))}
        else if((y[i]==4)) {m[i]<-log((1-pnorm(cut4.OP-ZG.OP))/((pnorm(XB.ZIOP))*((1-pnorm(cut4.ZIOP-ZG.ZIOP) ))))}     
	}
diff.m.sqrd<-(m-mean(m))^2
sum.d.m.s<-sum(diff.m.sqrd)
Vuong.OP.ZIOP<-((sqrt(n))*(1/n)*sum(m))/sqrt((1/n)*(sum.d.m.s))
Vuong.OP.ZIOP


#Vuong test comparing Main ZIOPC to Non-Truncated OP
m<-matrix(0,nrow=n,ncol=1)
  for(i in 1:n){					              					      
	ZG.OP<-(1.18620*z[i,1]+-0.01618*z[i,2]+-0.32590*z[i,3]+0.01053*z[i,4]+-0.14820*z[i,5])
	cut1.OP<-2.72527
	cut2.OP<-2.72527+exp(-4.28509)
	cut3.OP<-2.72527+exp(-4.28509)+exp(-2.46453)
	cut4.OP<-2.72527+exp(-4.28509)+exp(-2.46453)+exp(-0.65400)
	XB.ZIOPC<-(-2.76473*x[i,1]+3.59316*x[i,2]+0.29423*x[i,3]+-0.15366*x[i,4]+1.04496*x[i,5])
	ZG.ZIOPC<-(-0.78682*z[i,1]+-0.20957*z[i,2]+-0.44798*z[i,3]+0.09295*z[i,4]+0.06872*z[i,5])
	cut1.ZIOPC<-0.74757
	cut2.ZIOPC<-0.74757+exp(-3.77804)
	cut3.ZIOPC<-0.74757+exp(-3.77804)+exp(-1.97628)
	cut4.ZIOPC<-0.74757+exp(-3.77804)+exp(-1.97628)+exp(-0.20074)
	rho<-0.10578        
	negsigma<-matrix(c(1,-rho,-rho,1), nrow=2, ncol=2)
	sigma<-matrix(c(1,rho,rho,1), nrow=2, ncol=2)
        if((y[i]==0)){m[i]<-log(pnorm(cut1.OP-ZG.OP)/((1-pnorm(XB.ZIOPC))+pmvnorm(lower=c(-Inf,-Inf), upper=c(XB.ZIOPC,(cut1.ZIOPC-ZG.ZIOPC)),mean=c(0,0),sigma=negsigma)))}
        else if((y[i]==1)){m[i]<-log((pnorm(cut2.OP-ZG.OP)-pnorm(cut1.OP-ZG.OP))/(pmvnorm(lower=c(-Inf,-Inf), upper=c(XB.ZIOPC,(cut2.ZIOPC-ZG.ZIOPC)),mean=c(0,0),sigma=negsigma)-pmvnorm(lower=c(-Inf,-Inf), upper=c(XB.ZIOPC,(cut1.ZIOPC-ZG.ZIOPC)),mean=c(0,0),sigma=negsigma)))}
        else if((y[i]==2)){m[i]<-log((pnorm(cut3.OP-ZG.OP)-pnorm(cut2.OP-ZG.OP))/(pmvnorm(lower=c(-Inf,-Inf), upper=c(XB.ZIOPC,(cut3.ZIOPC-ZG.ZIOPC)),mean=c(0,0),sigma=negsigma)-pmvnorm(lower=c(-Inf,-Inf), upper=c(XB.ZIOPC,(cut2.ZIOPC-ZG.ZIOPC)),mean=c(0,0),sigma=negsigma)))}
        else if((y[i]==3)){m[i]<-log((pnorm(cut4.OP-ZG.OP)-pnorm(cut3.OP-ZG.OP))/(pmvnorm(lower=c(-Inf,-Inf), upper=c(XB.ZIOPC,(cut4.ZIOPC-ZG.ZIOPC)),mean=c(0,0),sigma=negsigma)-pmvnorm(lower=c(-Inf,-Inf), upper=c(XB.ZIOPC,(cut3.ZIOPC-ZG.ZIOPC)),mean=c(0,0),sigma=negsigma)))}
        else if((y[i]==4)) {m[i]<-log((1-pnorm(cut4.OP-ZG.OP))/(pmvnorm(lower=c(-Inf,-Inf), upper=c(XB.ZIOPC,(ZG.ZIOPC-cut4.ZIOPC)),mean=c(0,0),sigma=sigma)))}     
	}
diff.m.sqrd<-(m-mean(m))^2
sum.d.m.s<-sum(diff.m.sqrd)
Vuong.OP.ZIOPC<-((sqrt(n))*(1/n)*sum(m))/sqrt((1/n)*(sum.d.m.s))
Vuong.OP.ZIOPC

   

#Next, compare robustness ZiOP(C) models to non-truncated OP model

#set the parameters
  y<-newdata.onset[,13]		# hostility
  x1<-newdata.onset[,6]		# contig
  x2<-newdata.onset[,7]		# alliance
  x3<-newdata.onset[,8]		# dem
  x4<-newdata.onset[,9]		# power		                                             
  x5<-newdata.onset[,11]	# maturity one
  x6<-newdata.onset[,12]	# maturity none
  x<-cbind(1,x1,x2,x3,x4,x5,x6) # inflation stage covariates                                                   	                                             					     	      
  z<-cbind(x1,x2,x3,x5,x6)      # outcome stage covariates 
  n<-nrow(newdata.onset)	# number of observations
     
#Vuong test comparing Robustness ZIOP to Non-Truncated OP
m<-matrix(0,nrow=n,ncol=1)
  for(i in 1:n){				              					      
	ZG.OP<-(1.18620*z[i,1]+-0.01618*z[i,2]+-0.32590*z[i,3]+0.01053*z[i,4]+-0.14820*z[i,5])
	cut1.OP<-2.72527
	cut2.OP<-2.72527+exp(-4.28509)
	cut3.OP<-2.72527+exp(-4.28509)+exp(-2.46453)
	cut4.OP<-2.72527+exp(-4.28509)+exp(-2.46453)+exp(-0.65400)
	XB.ZIOP<-(-2.71893*x[i,1]+3.67429*x[i,2]+0.27915*x[i,3]+-0.17241*x[i,4]+0.98522*x[i,5]+0.05455*x[i,6]+-0.08415*x[i,7])
	ZG.ZIOP<-(-1.03160*z[i,1]+  -0.21248*z[i,2]  +-0.43051*z[i,3]+0.05248*z[i,4]+0.10722*z[i,5])
	cut1.ZIOP<-0.49784
	cut2.ZIOP<-0.49784+exp(-3.78138)
	cut3.ZIOP<-0.49784+exp(-3.78138)+exp(-1.97942)
	cut4.ZIOP<-0.49784+exp(-3.78138)+exp(-1.97942)+exp(-0.20222)
        if((y[i]==0)){m[i]<-log(pnorm(cut1.OP-ZG.OP)/((1-pnorm(XB.ZIOP))+(pnorm(XB.ZIOP))*(pnorm(cut1.ZIOP-ZG.ZIOP))))}
        else if((y[i]==1)){m[i]<-log((pnorm(cut2.OP-ZG.OP)-pnorm(cut1.OP-ZG.OP))/((pnorm(XB.ZIOP))*((pnorm(cut2.ZIOP-ZG.ZIOP)) - (pnorm(cut1.ZIOP-ZG.ZIOP)))))}
        else if((y[i]==2)){m[i]<-log((pnorm(cut3.OP-ZG.OP)-pnorm(cut2.OP-ZG.OP))/((pnorm(XB.ZIOP))*((pnorm(cut3.ZIOP-ZG.ZIOP)) - (pnorm(cut2.ZIOP-ZG.ZIOP)))))}
	else if((y[i]==3)){m[i]<-log((pnorm(cut4.OP-ZG.OP)-pnorm(cut3.OP-ZG.OP))/((pnorm(XB.ZIOP))*((pnorm(cut4.ZIOP-ZG.ZIOP)) - (pnorm(cut3.ZIOP-ZG.ZIOP)))))}
        else if((y[i]==4)) {m[i]<-log((1-pnorm(cut4.OP-ZG.OP))/((pnorm(XB.ZIOP))*((1-pnorm(cut4.ZIOP-ZG.ZIOP)))))}     
	}
diff.m.sqrd<-(m-mean(m))^2
sum.d.m.s<-sum(diff.m.sqrd)
Vuong.OP.ZIOP<-((sqrt(n))*(1/n)*sum(m))/sqrt((1/n)*(sum.d.m.s))
Vuong.OP.ZIOP




#Vuong test comparing Robustness ZIOPC to Non-Truncated OP
m<-matrix(0,nrow=n,ncol=1)
  for(i in 1:n){				              					      
	ZG.OP<-(1.18620*z[i,1]+-0.01618*z[i,2]+-0.32590*z[i,3]+0.01053*z[i,4]+-0.14820*z[i,5])
	cut1.OP<-2.72527
	cut2.OP<-2.72527+exp(-4.28509)
	cut3.OP<-2.72527+exp(-4.28509)+exp(-2.46453)
	cut4.OP<-2.72527+exp(-4.28509)+exp(-2.46453)+exp(-0.65400)
	XB.ZIOPC<-(-2.75717*x[i,1]+3.46711*x[i,2]+0.29482*x[i,3]+-0.17726*x[i,4]+1.03809*x[i,5]+0.05047*x[i,6]+-0.09596*x[i,7])
	ZG.ZIOPC<-(-0.70669*z[i,1]+-0.20976*z[i,2]+-0.43847*z[i,3]+0.05755*z[i,4]+0.11162*z[i,5])
	cut1.ZIOPC<-0.83365
	cut2.ZIOPC<-0.83365+exp(-3.76988)
	cut3.ZIOPC<-0.83365+exp(-3.76988)+exp(-1.96863)
	cut4.ZIOPC<-0.83365+exp(-3.76988)+exp(-1.96863)+exp(-0.19655)
	rho<-0.15128       
	negsigma<-matrix(c(1,-rho,-rho,1), nrow=2, ncol=2)
	sigma<-matrix(c(1,rho,rho,1), nrow=2, ncol=2)
        if((y[i]==0)){m[i]<-log(pnorm(cut1.OP-ZG.OP)/((1-pnorm(XB.ZIOPC))+pmvnorm(lower=c(-Inf,-Inf), upper=c(XB.ZIOPC,(cut1.ZIOPC-ZG.ZIOPC)),mean=c(0,0),sigma=negsigma)))}
        else if((y[i]==1)){m[i]<-log((pnorm(cut2.OP-ZG.OP)-pnorm(cut1.OP-ZG.OP))/(pmvnorm(lower=c(-Inf,-Inf), upper=c(XB.ZIOPC,(cut2.ZIOPC-ZG.ZIOPC)),mean=c(0,0),sigma=negsigma)-pmvnorm(lower=c(-Inf,-Inf), upper=c(XB.ZIOPC,(cut1.ZIOPC-ZG.ZIOPC)),mean=c(0,0),sigma=negsigma)))}
        else if((y[i]==2)){m[i]<-log((pnorm(cut3.OP-ZG.OP)-pnorm(cut2.OP-ZG.OP))/(pmvnorm(lower=c(-Inf,-Inf), upper=c(XB.ZIOPC,(cut3.ZIOPC-ZG.ZIOPC)),mean=c(0,0),sigma=negsigma)-pmvnorm(lower=c(-Inf,-Inf), upper=c(XB.ZIOPC,(cut2.ZIOPC-ZG.ZIOPC)),mean=c(0,0),sigma=negsigma)))}
        else if((y[i]==3)){m[i]<-log((pnorm(cut4.OP-ZG.OP)-pnorm(cut3.OP-ZG.OP))/(pmvnorm(lower=c(-Inf,-Inf), upper=c(XB.ZIOPC,(cut4.ZIOPC-ZG.ZIOPC)),mean=c(0,0),sigma=negsigma)-pmvnorm(lower=c(-Inf,-Inf), upper=c(XB.ZIOPC,(cut3.ZIOPC-ZG.ZIOPC)),mean=c(0,0),sigma=negsigma)))}
        else if((y[i]==4)) {m[i]<-log((1-pnorm(cut4.OP-ZG.OP))/(pmvnorm(lower=c(-Inf,-Inf), upper=c(XB.ZIOPC,(ZG.ZIOPC-cut4.ZIOPC)),mean=c(0,0),sigma=sigma)))}     
	}
diff.m.sqrd<-(m-mean(m))^2
sum.d.m.s<-sum(diff.m.sqrd)
Vuong.OP.ZIOPC<-((sqrt(n))*(1/n)*sum(m))/sqrt((1/n)*(sum.d.m.s))
Vuong.OP.ZIOPC



############################################################################
######################Likelihood Ratio Tests################################
############################################################################

#likelihood ratio test comparing the two ZiOP models
l.ratio<--2*((-11901)+11892)
l.ratio
pchisq(l.ratio, 2)


#likelihood ratio test comparing the two ZiOPC models
l.ratio<--2*((-11901)+11891)
l.ratio
pchisq(l.ratio, 2)









